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ABSTRACT 

I use cosmological simulations that incorporate a physically motivated 
approximation to three-dimensional radiative transfer that recovers correct 
asymptotic ionization front propagation speeds for some cosmologically relevant 
density distributions to investigate the process of the reionization of the universe 
by ionizing radiation from proto-galaxies. Reionization proceeds in three stages and 
occupies a large redshift range from z ~ 15 until z ~ 5. During the first, "pre-overlap" 
stage, H II regions gradually expand into the low density IGM, leaving behind neutral 
high density protrusions. During the second, "overlap" stage, that occurs in about 
10% of the Hubble time, Hii regions merge and the ionizing background rises by a 
large factor. During the third, "post-overlap" stage, remaining high density regions are 
being gradually ionized as the required ionizing photons are being produced. 

Residual fluctuations in the ionizing background reach significant (more than 10%) 
levels for the Lyman-alpha forest absorption systems with column densities above 
lO^'' - 10^^ cm~2 at z = 3 to 4. 



Subject headings: cosmology: theory - cosmology: large-scale structure of universe 
galaxies: formation - galaxies: intergalactic medium 



1. Introduction 



Existing ground based observations of the CMB on sub-degree angular scales suggest that 
the gas content of the universe was mostly neutral since recombination at z ^ 1000 until about 
z ~ 100 (Bond & Jaffe 1998; Griffiths, Barbosa, & Liddle 1998), because earlier reionization would 
have brought the last scattering surface to lower redshift, smoothing the intrinsic CMB anisotropy. 
At the same time we know that the universe is highly ionized since z ~ 5, from observations of the 
spectra of quasars with the highest redshifts (Giallongo et al. 1994; Williger et al. 1994; Songaila 
et al. 1999). This change of the ionization state of the universe from neutral to highly ionized is 
called reionization. 

Recent years witnessed a surge in research on reionization along two separate directions. 
Semi-analytical methods attempt to describe the general features of the evolution of the 
intergalactic medium (IGM) based on simple, ad hoc assumptions about the star formation history 



Gnedin: Cosmological Reionization 



2 



and the density distribution in the IGM. Most of the previous works utihzed a simple clumping 
factor to account for the inhomogeneity of the IGM (Giroux & Shapiro 1996; Tegmark et al. 1997; 
Madau, Meiksin, & Rees 1997; Ciardi & Ferrara 1997; Haiman & Loeb 1997, 1998; Madau, Haardt, 
k Rees 1999; Valageas & Silk 1999; Chiu & Ostriker 1999). A recent work by Miralda-Escude, 
Haehnelt, &; Rees (1999) attempted to advance further by introducing the distribution of the 
ionized fraction as a function of gas density which, besides everything else, provides a specific 
model for various clumping factors. Despite using a very simple (perhaps even oversimplified) 
ansatz for the ionized fraction - density distribution, their main conclusions agree remarkably well 
with the results of the simulations presented in this paper, as I will show below. 

The main advantage of the semi-analytical approaches is their simplicity and the ability to 
emphasize the key physical processes. Their main problem is a significant oversimplification and 
inability to follow the complex dynamical interactions between the dark matter, gas, stars, and 
the spatially inhomogeneous and time-variable radiation field. Semi-analytical models also miss 
the effects of correlation between the sources of ionizing radiation. 

The latter deficiency was partly overcome recently by Ciardi et al. (1999) by combining 
a semi-analytical model with the A?^-body simulation, but the lack of dynamics will elude 
semi-analytical models forever. 

Cosmological numerical simulations offer a totally different approach to modeling reionization 
(Ostriker &; Gnedin 1996; Gnedin & Ostriker 1997). The main advantage of simulations 
over the semi-analytical models is that numerical simulations fully account for the dynamical 
evolution of the matter contents of the universe, thus avoiding the main over-simplification of 
the semi-analytical models. Their inherited limitation is a limited dynamic range, which for 
existing numerical codes implies an impossibility to achieve numerical convergence and therefore 
quantitatively accurate results. However, semi-analytical models are also imprecise since they 
adopt ad hoc assumptions and ansatzes, and thus we must accept the fact that we are not able 
to model reionization quantitatively accurate at this moment, and need to limit ourselves with 
understanding the key qualitative features of the physical processes that take place in the universe 
during that epochs. 

Until recently, numerical simulations failed to include the three-dimensional radiative transfer 
in all its complexity, and therefore were not sufficient to properly model the dynamics of the 
ionizing radiation in the universe. 

The understanding of the crucial role that the full 3D radiative transfer plays in the physics 
of reionization prompted the first attempts to develop numerical techniques to incorporate the 
radiative transfer into numerical simulations (Abel, Norman, &; Madau 1999). While Abel et al. 
(1999) were able to develop an algorithm for an exact implementation of the full 3D radiative 
transfer into cosmological simulations, their algorithm is appropriate for following the radiation 
field in the case when there are only a few sources within the volume with the size of the mean 
free path of an ionizing photon. When the number of sources within this volume increases, as is 
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the case after all Hii regions overlap, their method becomes prohibitively expensive. 

In this paper I adopt a somewhat different approach. Rather than trying to model the 
radiative transfer exactly, I develop an approximation, which is sufficiently fast computationally. I 
validate the adopted approximation on a series of spherically symmetric test cases, for which the 
full radiative transfer can be solved exactly at a modest computational expense. 

The physical basis for the adopted approximation, which I call a "Local Optical Depth" 
(LOD) approximation, is simple. In the optically thin regime it is straightforward to compute the 
spatially variable radiation field, as the problem reduces to simply collecting potential laws 
from all the sources in the simulation volume, which can be done fast with the P^M algorithm. 
Thus, the key difficulty in following the full radiative transfer is to account for the optical 
depth, which can always be presented as a product of the cross-section, the gas density, and a 
characteristic scale. The LOD approximation adopts an ansatz for this characteristic scale as the 
scale over which the gas density changes significantly. 

Having developed the means to follow approximately the 3D spatially inhomogeneous and 
time-dependent radiative transfer, I apply them to modeling the process of reionization of 
the universe in a specific scenario when sources of ionizing radiation are stars (grouped into 
proto-galaxies) . There are two reasons why numerical simulations need to be restricted to this 
specific case at the moment. 

The first reason is "scientific" , as the mounting evidence suggests that quasars are unable to 
reionize the universe at 2 > 5 (Madau et al. 1999), and thus the universe is reionized by stars. The 
second reason is "pragmatic". A quasar radiation field may affect the properties of the IGM over 
a region several tens of megaparsecs in size. Thus, a numerical simulation capable of modeling the 
reionization by quasars would have to have a computational region of at least 100 Mpc in size. At 
the same time, in order to resolve the formation of the first structures and to follow the gas density 
with sufficient precision, it will need spatial resolution below 1 kpc and mass resolution better 
than about 10^ solar mass in baryons. Such a dynamic range (10^) is beyond the capabilities of 
the existing cosmological hydrodynamic codes.[| 

Thus, only a numerical simulation with the box size of several Mpc is practical at the 
moment. Fortunately, this box size is sufficient to model the reionization by stellar sources (on a 
semi-qualitative level), as will be shown below. 

This paper is organized in a conventional way. In §2 I describe the simulations, §3 is devoted 
to the results, and §4 concludes with the discussion. Appendix describes the Local Optical Depth 
approximation and the tests of the method. 



^With the exception of the Adaptive Mesh Refinement technique. 
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2. Simulations 

Simulations reported in this paper were performed with the "Softened Lagrangian 
Hydrodynamics" (SLH-P^M) code (Gnedin 1995, 1996; Gnedin k Bertschinger 1996). The 
following physical ingredients are currently included in the code: 

Dark matter is followed using the adaptive P^M algorithm. 

Gas dynamics is followed on a quasi-Lagrangian deformable mesh using the SLH algorithm. 

Star formation is included using the Schmidt law in resolution elements that sink below the 
numerical resolution of the code. 

Atomic physics of hydrogen and helium plasma is followed exactly using a two-level implicit 
scheme. 

Molecular hydrogen formation and destruction is followed exactly (including the radiative 
transfer effects) in the limit when the fraction of hydrogen in the molecular form is small 
and the self-shielding of H2 is unimportant (in the latter case the approximate method of 
following the radiative transfer becomes exact). This is always the case in the simulation 
presented in this paper because the numerical resolution is not sufficient to resolve the 
formation of molecular clouds. 

Radiative transfer is treated self-consistently, albeit only approximately, in a 3D spatially- 
inhomogeneous and time-dependent manner. This is a new addition to the SLH-P^M code. 
The specific implementation of the radiative transfer is described in the Appendix. 

Hereafter I adopt the following values of cosmological parameters in the CDM-I-A model: 

no = 0.3, = 0.7, h = 0.7, = 0.04. 

This cosmological model is in a reasonable agreement with all available observational data and 
has a desirable feature that all the parameters are specified to only one decimal place. For the 
simulations presented in this paper the particular choice of cosmological parameters matters rather 
little as all reasonable models look similar at high redshift. 

There are two free parameters that control the behavior of the simulations. The first 
parameter is the efficiency of star formation e*, that appears in the Schmidt law: 

dt ~'*t,' 

where and pg are the stellar and the gas density respectively, and is the maximum of 
the dynamical and cooling time. Equation (|l]) is only applied in resolution elements that are 
determined to be beyond the resolution limit of the simulation. In all resolved resolution elements 
star formation is not allowed. 
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The second important parameter measures the amount of ionizing radiation each "star" emits: 

duuv 2 dp* 

^ = ^^^^ ^' 

where uuv is the production of the energy density in ionizing radiation per unit time. This 
parameter apparently depends on the initial mass function (IMF). For the Salpeter IMF the 
"escape fraction" of the ionizing photons, i.e. the fraction of ionizing photons that escape from the 
immediate vicinity of a star, is roughly 

/esc ~ 1.4 X W^UV (3) 

(Madau et al. 1999). Thus, an efficiency of 4 x 10^^ corresponds to roughly 60% escape fraction. 

However, what matters most for the evolution of the IGM is only the total number of ionizing 
photons emitted, i.e. the product of e* and euv- Therefore, I fix = 0.05 as my fiducial value, 
which is similar to the values found in the Milky Way. The latter by itself is not a justification for 
this particular choice, since there may be very little in common between the star formation in the 
Milky Way and at high redshift. However, this choice leads to about 3 to 5 percent of all baryons 
being converted to stars by z = 4, which translates to about 10 to 20 percent of all stars being 
formed by 2: = 4, a reasonable estimate given the current data on the star formation history of the 
universe (Madau 1999, Steidel et al. 1999, Renzini 1999). In addition, the star formation rate at 
z = 4 in the simulation turns out to be similar to the observational value (see Fig. ^ below). 

It is important to emphasize here that an assumed value of the radiation efficiency, or for that 
matter of the escape fraction, is resolution dependent, since it measures the fraction of energy that 
escapes from simulated stars into the simulated gas. Since any simulation has a finite resolution, 
this means that the escape fraction in a simulation measures the fraction of radiation that 
"escapes from the stellar surface to the resolution scale of the simulation", and mathematically 
precise meaning of this phrase cannot be formulated, as it may depend on the behavior of a 
numerical scheme at its resolution limit. It may even be a bad approximation to assume ejjv to be 
constant during the simulation, since simulations presented in this paper have a fixed resolution 
in comoving coordinates, which means that in physical units the resolution length increases with 
time. However, since the process of reionization takes a relatively short period of time, this effect 
is likely to be insignificant. At any rate, as I argue below, the simulations presented in this paper 
are only reliable on a semi-qualitative level (within a factor of two or so), and cannot be relied 
upon to provide an accurate quantitative answer. 

Table || lists all the simulations presented in this paper. The last column in the table is the 
final redshift at which a given simulation is stopped. The first four simulations whose labels begin 
with N64 contain 64^ dark matter particles, 64^ baryonic mesh, and a number of stellar particles 
which form during the simulation. These runs have a dynamic range of 1300. All 64^ runs have a 
computational box size of 2h~^ comoving megaparsecs. 

The 64^ runs serve as tests and designed to investigate the parameter dependence and the 
convergence of the simulations. They are not used for producing scientific results. The first three 
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Table 1. Numerical Parameters 



Run 


N 


Box size 


Baryonic mass res. 


Spatial res. 


e* 








N64_L2_A 


64^ 


2/1-1 ]YIpc 


lO^-'^ M0 


1.5/1-1 kpc 


0.05 


4 X 10- 


5 


4 


N64_L2_B 


64^ 


2/i-i Mpc 




1.5/1-1 kpc 


0.05 


1.2 X 10 


-4 


4 


N64_L2_C 


64^ 






1.5/1-1 kpc 


0.05 


1.3 X 10 


-5 


4 


N64_L2_D 


64^ 


2/1-1 Mpc 




1.5/1-1 kpc 


0.10 


4 X 10- 


5 


4 


N128_L4_A 


128^ 


4/i-i Mpc 




l.O/i-i kpc 


0.05 


4 X 10- 


5 


4 


N128_L8_A 


128^ 


8/i-i Mpc 




2.0/1-1 kpc 


0.05 


4 X 10" 


5 


6.5 


N128_L2_A 


128^ 


2/i-i Mpc 


W^-^ Mq 


0.5/1-1 j^pp 


0.025 


4 X 10- 


5 


6.5 



of them differ only by the value of the radiation efficiency, and the last run, N64_L2_D, has a twice 
higher star formation rate. One may note that the escape fraction in run N64_L2_B is formally 
170% if equation (^) is used. This should not be considered unphysical since this run is only used 
as a test, and besides that a 70% increase in the fraction of high mass stars will compensate for 
that. Anyway, as I show below, this value of the radiation efficiency does not fit the observations. 

The last three simulations listed in Table [l| whose labels begin with N128 incorporate 128'^ 
dark matter particles, the same number of baryonic cells, and about twice more stellar particles, 
which keep forming during the simulation. The dynamic range of these three simulations is 4000, 
and other parameters are listed in the table. Run N128_L4_A is the production run of this paper, 
and two other 128^ runs are also used to investigate numerical convergence, and therefore are not 
continued until z = 4. 

Figure || illustrates in a graphical form the mass and spatial scales simulated in this paper. 
The production run N128_L4_A has just enough mass resolution to resolve the filtering scale, the 
minimum scale on which baryons cluster in linear theory.^ This is clearly only marginally enough, 
as one would prefer to resolve the filtering scale by a factor of 10 or more to have a safe margin. 
On the other hand the production run N128_L4_A does not miss a large amount of small scale 
power. Thus, even from this consideration we should expect that the simulations presented in 
this paper are not quantitatively accurate, but should be sufficient to give a general qualitative 
picture, since they include all the relevant physics and account for most of the initial power. 

Figure |2| now demonstrates the level of convergence of the simulations presented in this paper 
and the method of choosing the parameter ejjv- Three thin lines: dotted, short-dashed, and 



^As has been shown in Gnedin & Hui (1998) at z = 10 the mass corresponding to this scale is about 11 times 
larger than the Jeans mass at this redshift. 
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Fig. 1. — Mass and spatial resolution of the simulations presented in Table |l|. Left y axis labels comoving 
scales and right y-axis labels physical scales at z = 10. The thin solid box shows all four 64^ runs (slightly 
offset for clarity), the bold solid box shows run N128_L4_A, and two bold dashed boxes show runs N128_L2_A 
and N128i8_A. The shaded regions marks the mass scales on which baryons do not cluster in linear theory 
due to finite pressure (the so-called filtering scale; Gnedin & Hui 1998). 

long-dashed show runs N64_L2_C, N64_L2_A, and N64_L2_B respectively, which only differ by the 
value of ejjv One can see that at z = 4 the value of the ionizing intensity J21 is 0.06, 0.4, and 
2.5, respectively, changing by a factor of about 6.5 whenever the radiation efficiency changes by a 
factor of 3. Since the observational values for J21 at z = 4 range between 0.2 and 0.4 (Lu et al. 
1996), I choose ejjv = 4 x 10^^ (used in run N64_L2_A) as a reasonable value for this parameter. 
Since my simulations cannot be used for accurate quantitative predictions, such a rough agreement 
with the observations is sufficient. 

A similar exercise with 128^ simulations is however not feasible, as it would require an 
excessive amount of computer time to run three 128^ simulations up to z = 4 with three different 
values of ejjv However, as can be seen from the top panel of Fig. |2|, the production run simulation 
N128_L4_A with assumed value of euv = 4 x 10~^ also gives J21 = 0.4 at z = 4, and thus a guess 
based on 64^ simulations turns out to be a good one, whereas a small run N64_L2_D, which has 
twice higher star formation rate than N64_L2_A, agrees better with the production run before and 
during the epoch of overlap (indicated by a sharp rise in J2i). 

The two other 128^ runs were normalized so as to agree with the production run N128_L4_A 
at z ~ 10 in the star formation rate and the ionizing intensity J21 (the small box run N128_L2_A 
actually has somewhat lower star formation rate). Nevertheless, they disagree with it at the epoch 
of overlap, and this difference indicates the lack of numerical convergence. 
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Naeamme et al 1999 _ 




Steidel et al 1999 



Fig. 2b 



Fig. 2. — Tlie evolution of the spatially averaged ionizing intensity J21 (measured in conventional units 
of 10^^^ erg/ cm^/ sec/ Hz/ rad) (a) and the star formation rate (&) as a function of redshift for seven runs 
listed in Table |^: N64_L2_A {thin short-dashed line), N64_L2_B {thin long-dashed line), N64X2_C {thin 
dotted line), N64_L2_D {thin solid line), N128_L4_A {bold solid line), N128_L8_A {bold long-dashed line), 
N128i2_A {bold dotted line). The data points for the star formation rate at z « 4.2 are from Steidel et 
al. (1999; filled circle) and the same data corrected for dust extinction by a different method by Nagamine, 
Cen, & Ostriker (1999; empty circle). 

Let me first focus on an 8h~^ Mpc run N128_L8_A. It agrees quite well with the production 
(4/i~^ Mpc) run before the overlap, but predicts a somewhat lower redshift of overlap. It however 
agrees better with the smaller N64_L2_A run, because its spatial resolution is closer to the small 
run than to the production run. Thus, the redshift of overlap is delayed if the spatial resolution of 
a simulation is not sufficient. Is the resolution of the production run sufficient in this respect? We 
can use comparison between the production run N128_L4_A and a 2h~^ Mpc run N128_L2_A to 
investigate the eff'ect of spatial resolution. It is clear that the production run still lacks resolution 
to resolve the very first star formation at z ~ 20. This is not surprising, since it barely resolves 
the baryon filtering scale before the reheating, and thus misses about 7% of the total small scale 
power. At z ~ 10 the two runs more-or-less agree, but the smaller run has a much earlier time 
of overlap. Also, it gives J21 at z = 6.5 of about 0.6, already higher than the final value of the 
production run of about 0.4. Thus, run N128_L2_A can be considered as giving the upper limit 
for the redshift of overlap for this model (which is 2; = 8.5), and this upper limit is likely to be 
too high, because run N128_L2_A is going to have a too high a value for J21 at lower redshifts. 
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This is of course not surprising, since run N128_L2_A has the same box size as N64_L2_A, but 
much higher spatial and mass resolution, so it is going to have more star formation than a smaller 
run, and thus an earlier reionization. In comparison to the production run N128_L4_A, it misses 
important large scale waves which contribute toward absorption. It also suffers from the cosmic 
variance problem, since at z = 4 the rms density fluctuation in a 2h~^ Mpc box is 0.41, whereas 
for a 4/i-i Mpc it is 0.27, and for a 8/i"^ Mpc it is 0.16. 

The comparison above does not however constitute the complete test of numerical convergence, 
as three different resolutions are required for such a test. Thus, a 256^ run is needed to assign a 
meaningful value for the numerical error of my simulations, but such a run is currently beyond 
the available means. I must therefore acknowledge that the simulations presented in this paper 
can only be considered reliable on a semi-qualitative basis, within a factor of two or so. One can 
use the difference between the bold solid line and thin solid and dashed lines as an estimate of 
the numerical error, but without a larger simulation it is impossible to assign any meaningful 
confidence level to this error. 

3. Results 

3.1. Reionization at a glance 

We can now look at a general description of the reionization process, as refiected in the 
simulations. In order to visualize the process, I show in Figure ^ a thin {15h~^ comoving 
kiloparsecs deep) slice through the computational box taken at a random place. The three panels 
show in a logarithmic stretch the neutral hydrogen (upper left panel), the gas density (lower left 
panel), and the gas temperature (lower right panel). The upper right panel shows the evolution of 
the ionizing intensity (logarithmically scaled) versus the redshift (the same as the solid bold line 
in Fig. ^). The scale bars next to respective panels show the correspondence between the color 
and the numeric scale. The epoch as measured by the scale factor a = 1/(1 + z) is labeled at the 
center of each plot.^ 

I can now describe the general features of reionization. The reionization starts (Fig. ^) 
with ionization fronts propagating from proto-galaxies located in high density regions into the 
voids, leaving the high density outskirts of an object still neutral, because at high density the 
recombination time is very short, and there are not enough photons to ionize the high density 
regions. As an ionization front moves on with the Hii region forming behind it, it leaves behind 
high density regions which require many more photons to get ionized than is available at that 
moment (Fig. ^-d). This stage of the reionization process can be called "pre-overlap" , and it 
extends over a considerable range of redshifts Az ~ 5 around z ~ 10. During this time the high 



^An MPEG movie of the simulation is available at 



http : //casa. Colorado . edu/'^gnedin/GALLERY/rei.p .html. 
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Fig. 3a 



Fig. 3b 





Fig. 3c 



Fig. 3d 



Fig. 3. — A thin slice througli tlie simulation volume at eight different epochs: (a) z — ff.5, (6) z = 9, 
(c) z — 7.7, (d) z — 7, (e) z — 6.7, (/) z — 6.1, (g) z — 5.7, (h) z — 4.9. Shown are logarithm of neutral 
hydrogen (upper-left), logarithm of gas density (lower-left) , logarithm of gas temperature (lower-right), and 
lg(J2i) as a function of redshift (upper-right). 

density regions around the source are slowly becoming ionized, whereas high density regions far 
from the source remain neutral. The ionizing intensity at this time remains low and is slowly 
increasing with time. However, since the radiation field is highly inhomogeneous at this time, the 
ionizing intensity J21, which is by definition a space-average, has only a formal meaning. 



By 2; ~ 7 the Hii regions start to overlap (Fig. |3|d), with the result that the number of sources 
shining at an average place in the universe increases, and the ionizing intensity starts to rise 
rapidly. The process of reionization enters its second stage, the "overlap" , which is quite rapid 
(Fig. 0d-e). As the ionizing intensity is rapidly increasing, the last remains of the neutral low 
density IGM are quickly eliminated, the mean free path increases by some two orders of magnitude 
(as explained below) over a Hubble time or so, and voids become highly ionized (neutral fraction 
of the order of 10~^). The high density regions at this moment are still neutral, as the number of 
ionizing photons available is not sufficient to ionize them. 
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Fig. 3g Fig. 3h 



After the overlap is complete, the universe is left with highly ionized both low density regions 
and some of the high density ones (which happened to lie close to the source, where the local 
value of the ionizing intensity is higher than the spatial average). High density regions far from 
any source remain neutral. This stage can be called "post-overlap", and is well described by the 
semi-analytic model of Miralda-Escude et al. (1999). As time goes on, and more and more ionized 
photons are emitted, the high density regions are gradually being eaten away (Fig. ^-g), and the 
spatially averaged ionizing intensity J21 continues to rise slowly until it flattens out at z ~ 5 (Fig. 
^g). The latter effect is however likely to be an artifact of a finite simulation box, as by this time 
the mean free path exceeds the box size by about a factor of 10 (as is shown below), and the 
simulation runs out of high enough density peaks, which would be present in the real universe (or 
a larger box simulation), and fails to reproduce the reality even on a semi-qualitative level. 

Therefore, the recent argument of whether reionization is "fast" or "slow" (see, for example, 
Miralda-Escude et al. 1999 versus Gnedin & Ostriker 1997 or Madau et al. 1999) is, in large part, 
a question of terminology. If one considers the whole process of reionization, which consists of 
"pre-overlap" , "overlap", and "post-overlap", one inevitably concludes that reionization is slow 
(i.e. taking place over a Hubble time or longer). However, if one looks only at the process of 
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1+z 

Fig. 4a 



T3 




1+Z 



Fig. 4b 



Fig. 4. — The evolution of the mean free path (a) and its rate of change (&) for the production run N128i4_A 
{bold solid line) and two smaller runs N64_L2_D {thin solid line) and N64_L2_A {thin short-dashed line). Two 
thin dotted lines are explained in the text. 



"overlap" and labels it reionization, then one concludes that reionization is fast. 



3.2. Reionization in more detail 

Let me now focus on details of the process that is overviewed in the previous subsection. 
Figure ^ shows the evolution of the mean free path and its time derivative for the production 
run N128_L4_A and two smaller runs (the latter two to demonstrate the level of numerical 
convergence). The epoch of overlap is clearly distinguished by a sharp rise in the mean free path. 
The rate of this rise, some 100 times faster than the expansion of the universe, shows no trend of 
diminishing with the increase of the box size (the bold line versus thin ones), and thus is unlikely 
to be an artifact of a finite simulation box. 

Can this behavior be understood in simple terms? Let me consider a simple model for 
the universe: each object has a 1/r^ density profile around it until those profiles overlap with 
neighboring objects. A source, placed at the center of an object produces the photoionization rate 
that falls off as l/r"^ in the optically thin regime, 

r(r) = ro^. 
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Given the density profile with a core radius tq, 

/ \ '0 
n[r) = no— — 
Tq + 

the neutral hydrogen fraction at r ^ rg is constant, 

Rn(r) RriQ 

= = W = T7- 

The mean free path then is given by the following expression, 

where k is the photoionization rate averaged value of the absorption coefficient (I omit here 
frequency dependence for simplicity): 

(see Gnedin & Ostriker 1997, eqn. [A6]), and (F) is the photoionization rate spatially averaged 
over a sphere of radius R, 

(r) = A/;r(0,-^* = 3r„|. 

Using the expression for the density, for the mean free path I obtain: 

2 R 

\ = , 5 

where a is the ionization cross section. If, instead of I used a spatially averaged absorption 
coefficient (/c), which would be appropriate if the ionizing background was homogeneous, I would 
get a different result: 

c _ 1 R^ 
{k) 3crxono rl ' 

Reasonable values for uq and tq would be the virial density, no = 200nft = 4 x 10^^(1 + z)^ cm~^, 
where fii, is the average baryon density, and the virial radius, vq = 31(1 + z)^^ kpc for a 10^ 
solar mass object. Obviously, R has to be about the size of the Hii region around the source, 
which from Fig. ^ can be roughly approximated as i? = 10(1 + z)~^ Mpc at z 9. I also assume 
Xq ~ 10^^, which is justified below. 

The two dotted lines in Fig. |^ now show A (the lower line) and A (the upper line). Before the 
overlap, the radiation field is highly nonuniform, and the mean free path is short, as the majority 
of photons are absorbed within one Hll region. The lower dotted line (A from eqn. ^) provides in 
this regime a remarkably good fit to the simulation results despite the highly simplistic nature of 
the above derivation. Fig. ^p-d illustrate this effect. The Hii region located in the middle right 
of the slice contains the high density neutral filaments, which are slowly being eaten away by the 
ionization front, which moves much slower across high density filaments than over the low density 
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Fig. 5. — The joint mass-weighted distribution of the gas density and the local photoionization rate T in the 
production run N128_L4_A at four different redshifts. The straight bold line in the upper left panel marks 
a r = (r)(l + (5) law. 

void. The absorption in these high density regions shortens the mean free path compared to the 
size of the Hii region in the low density medium. This is consistent with the simplistic picture 
above which assumes a spherically symmetric density distribution, which in reality is an average 
of the high density filaments (with the density run shallower than 1/r^) and the low density voids 
(with the density run steeper than 1/r^). 

After the overlap, the ionizing radiation is more (albeit non perfectly) homogeneous, and 
the mean free path increases by almost two orders of magnitude and becomes comparable to the 
estimate A. Thus, the effect of the overlap being "fast" is entirely due to the change in the relative 
distribution of the ionizing radiation and the gas density, and thus, in essence, is a topological 
effect. 

This is further illustrated by Figure |5|, which shows on four panels the joint mass- weighted 
distribution of the gas density and the photoionization rate at four different redshift. At z = 9 the 
photoionization rate in the regions that are indeed photoionized (only those count towards the 
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Fig. 6. — The joint mass-weighted distribution of the gas density and neutral hydrogen fraction in the 
production run N128_L4_A at four different redshifts. The excessive sharpness of the features at (5 « is the 
defect of the radiative transfer approximation, in reahty there should be a more gradual transition from the 
underdense to the overdense regions. 

mean free path, since no photons are absorbed where there are no photons) is proportional to the 
gas density, with a large scatter reflecting the fact that different objects have different luminosities 
(in terms of the simplistic model above, Tq is different in different objects). As the Hii regions 
start to overlap (z = 7), regions with more or less uniform photoionization rate starts to appear 
at low densities. Finally, in the post-overlap stage, a large fraction of the total volume has a quite 
uniform photoionization rate, which still varies significantly in the high density regions {z = 6 
and z = 5). I will address the question of the inhomogeneity in the photoionization rate at low 
redshifts below. 

Figure ^ illustrates the process of reionization from yet another point of view. It presents 
the mass-weighted joint distribution of the gas density and neutral hydrogen fraction shown at 
the same four redshifts z = 9, 7, 6, and 5. In the pre-overlap stage the neutral hydrogen fraction 
is roughly independent of density in the photoionized regions, again with a large spread due to 
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different luminosities of different objects and with the median value of about xhi ~ 10~^, which 
is a justification for adopting this value in the simplistic estimate presented above. As time 
progresses and the Hi regions start to overlap, a regime xhi oc (1 + 5) starts to develop in the 
lower density gas, which is just a reflection of the fact that the photoionization rate is becoming 
independent of the density in those regions, whereas the high density gas still maintains the 
neutral fraction that is roughly independent of the density. After the overlap, most of the gas sits 
at the ionization equilibrium in a tight relation with the density, which becomes somewhat more 
spread in the highest density regions which correspond to the inner parts of individual objects. 
However, a considerable fraction of the high density gas remains in the neutral state, in accordance 
with the semi-analytical picture of Miralda-Escude et al. (1999). As the production of the ionizing 
photons continues, this gas is slowly being ionized, mostly the lower density regions first. 

One of the outcomes of the main assumption of this paper - that the sources of ionization are 
proto-galaxies rather than quasars - is that each halo harbors a source (except, perhaps, very low 
mass ones, that are not well resolved in my simulation), and therefore the halos are also ionized 
by local sources. The resolution of my simulation is however insufficient to resolve the interstellar 
medium, and proto-galaxies will contain giant molecular clouds which will be responsible for the 
damped Lyman-alpha absorption. 



3.3. What does it? 

The understanding of the process of reionization would not be complete without identifying 
which sources are indeed responsible for the bulk of ionizing photons. To this end, I have 
identified all bound objects in the production run N128_L4_A at z = 7 (the epoch of overlap) using 
Bertchinger & Gelb's (1991) DENMAX algorithm. Figure |^ shows the star formation rate versus 
the stellar mass for all objects. As can be seen, there exists a good correlation between the star 
formation rate and the mass of the object, 

M « ^ ^ 

* 1.6 X 108 Mq yr ' 

reflecting the fact that all objects found in the simulation form stars with the rapidly increasing 
rate (Fig. |). 

However, not all objects contribute equally toward the reionization of the low density IGM, 
but only those whose ionization fronts are able to leave their own halo. While the detailed analysis 
of the evolution of each individual ionization front is beyond the scope of this paper, a simple 
picture can be introduced that captures the main qualitative features of the simulation. 

If each individual object can be approximated as a homogeneous sphere with the radius equal 
the virial radius of the object, and the mean overdensity of 200, than an object will contribute 
towards the reionization of the low density IGM if the virial radius is smaller than the Stromgen 
radius of the ionization front. Figure ^b now shows the the Stromgen radius versus the virial 
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Fig. 7. — The star formation rate versus the stellar mass (a) and the virial radius versus the Stromgen 
radius (&) for all bound objects in the production run N128_L4_A. 

radius for all the bound resolved objects in the production run. One can see that the majority of 
objects have their Stromgen spheres well inside the virial radius, and only a handful of objects 
contribute to the ionizing flux that is capable of ionizing the low density IGM. 

This is not surprising, given what is seen in Fig. ^, and the fact that I have not been able to 
establish numerical convergence - if only a few objects in my simulation box are responsible for 
reionization, then the 4/i~^ Mpc box size is too small to be a representative volume of the universe, 
and a larger box size is needed to achieve the full numerical convergence. 

Figure ^ serves to illustrate this further. It shows the stellar mass function and the star 
formation rate function for the bound objects from my simulation, as well as the Schecter function 
fit to the stellar mass function (the fit deviates at low masses due to photoionization feedback). In 
addition, I plot with the dotted line the modified star formation rate function which includes only 
photons capable of ionizing the low density IGM, defined as 



Ml 



dn 

dMi' 



where 



Mi = max(0, {vs/r^f - 1)M* 
is proportional to the number of photons that escape into the IGM. As can be seen from Fig. it 
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Fig. 8. — The stellar mass function (solid bold line) and the star formation rate function (dashed line) at 
z = 7. The dotted line shows the modified star formation rate function which includes only photons capable 
of ionizing the low density IGM. The thin solid line is a Schecter function fit to the stellar mass function 
with a = -0.5 and = 3 x 10^ Mq. 

is objects with luminosities in excess of O.IL^ that are responsible for reionization, with and 
more luminous objects making the main contribution. 

3.4. Comparison with the semi-analytical models 

I have mentioned already that my simulations agree with the qualitative predictions of 
Miralda-Escude et al. (1999) quite well in the post-overlap stage. Figure ^ serves to illustrate the 
agreement even further. Fig. ^ shows a comparison between the mean volume emissivity e and 
the global recombination rate R, as introduced in equation (2) of Miralda-Escude et al. (1999). In 
accordance with their prediction, i? ~ e just after the overlap, and continues to be close to e and 
somewhat lower when other terms become important. 
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Fig. 9a Fig. 9b 

Fig. 9. — (a) Comparison of the mean volume emissivity e {solid line), measured in terms of the number of 
ionizing photons emitted for each atom in the universe per Hubble time, and the global recombination rate 
R (dashed line), measured in terms of the mean number of recombinations per Hubble time per baryon, as 
introduced in equation (2) of Miralda-Escude et al. (1999). (&) The time evolution of the quantity as 
defined in equation (1) of Miralda-Escude et al. (1999). 

Fig. |9|b shows the evolution of the quantity Aj defined in equation (1) of Miralda-Escude et 
al. (1999). Physically, Aj measures the cosmic density (in units of the mean density) above which 
the gas is mostly neutral, and below which the gas is mostly ionized. While Fig. ^ demonstrates 
that at high density the ionization state of the gas is not related to the density in any simple way, 
but also depends on whether a given fluid element is close or far from the nearest source, it is still 
possible to define Aj after the overlap as the limiting density above which 95% of the neutral gas 
lies. This quantity is plotted in Fig. ^b. While specific numerical values for Aj are significantly 
lower than those quoted in Miralda-Escude et al. (1999), because there is also the high density 
ionized gas responsible for a large number of recombinations, the general trend of Aj increasing 
with time is clearly observed. 

Figure |l^ shows the evolution of the mean number of ionized photons per baryon N^/i, on top 
and the neutral hydrogen fraction at the bottom for the production run and for the corresponding 
small run. One can see that by the time of overlap z = 7, about 10 ionizing photons were emitted 
per baryon, which may be considered a contradiction with Miralda-Escude et al. (1999). I must 
however reemphasize here that the total number of photons emitted depends on the radiation 
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Fig. 10. — Top panel: the mean number of ionizing photons emitted per baryon as a function of redshift 
for the production run N128_L4_A {bold line) and a smaUer run N64_L2_D {thin line line). Bottom panel: 
the mass- {dotted lines) and volume-averaged {solid lines) neutral hydrogen fractions for the production run 
N128_L4_A {bold lines) and a smaller run N64_L2_D {thin solid lines). 

efficiency euv, which in turn is highly resolution dependent. Therefore the quantity N.y/f, is also 
resolution dependent and has no physical meaning taken at a face value. The quantity used by 
Miralda-Escude et al. (1999) is however the number of ionizing photons per baryon escaped into 
the IGM from the immediate neighborhood of the ionizing sources ,0 and is therefore considerably 
smaller than the number directly taken from the simulation. The latter quantity cannot be much 
less than 1, and since it has to be considerably smaller than 10, the best estimate one can make is 
that it is of order of a few, again in agreement with Miralda-Escude et al. (1999). 

A large number of workers has studied reionization within the framework of the semi-analytical 
modeling using the clumping factor approach (see the references in the Introduction). Their 



^This quantity is physically intuitive, but is difficult to define in a mathematically rigorous fashion. 
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Fig. 11. — (a) The mean global recombination rate in units of the recombination rate of the fully ionized 
homogeneous universe (solid line), the density clumping factor (dotted line), and the ionized hydrogen 
clumping factor (dashed line) as a function of redshift. (b) The porosity of the IGM as a function of rcdshift 
for the approximation based on equation (23) of Madau et al. (1999) (dashed line) and the full calculation 
based on their equation (21) (solid line). 

main conclusion is that reionization is fast (because the recombination time is short), and that 
many ionizing photons per baryon are required to ionize the universe (again because of the same 
reason). It is therefore useful to compare this approach to the simulations and to Miralda-Escude 
et al. (1999). The latter comparison is especially instructive since there appears to be a complete 
disagreement between the two semi-analytical approaches: if Miralda-Escude et al. (1999) claim 
that reionization has to be gradual, the clumping factor approach predicts a fast reionization; if 
Miralda-Escude et al. (1999) show that just a few (perhaps as low as one) ionizing photons per 
baryon are needed for the overlap, the clumping factor approach requires that this number be 
large. 



I now show in Figure 11a the evolution of the volume averaged recombination rate 
(-R(T)nenHn) in units of the recombination rate for the fully ionized homogeneous universe with 

riH is the total number density of atomic and ionized hydrogen, 
and no = ran + nue, thus assuming that helium is only singly ionized). Also shown the density 
clumping factor 
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and the ionized hydrogen clumping factor 

{R{T)nenun) 



Ri{ne){nnu) ' 

which is simply related to the volume averaged recombination rate, 

{R(T)neniiu) = RineUnuCun, (7) 

where bar means volume average. 

There are few things to note. First, the recombination rate increases fast with time, and at 
the epoch of overlap already approaches 100 times that of the homogeneous universe. The density 
clumping factor is only slightly higher than the recombination rate, and the ionized hydrogen 
clumping factor begins its evolution with very high values, and approaches the recombination 
rate at later times, as can be expected from equation This behavior is quite different from 
what was predicted in Gnedin & Ostriker (1997), and the reason for the disagreement is quite 
clear: since the ionization fronts originate in the highest density regions, the first Hii regions are 
extremely overdense, which results in very large values for the clumping factor. However, since 
these first Hii regions occupy only a small fraction of the volume, the recombination rate is not 
large. 

Thus, the simulations confirm the predictions of the clumping factor approach: the clumping 
factor is indeed large at the time of overlap. To do a more detailed comparison, I picked up a 
Madau et al. (1999) paper as an example of semi-analytical modeling using this type of approach. 
There have been recently a series of detailed papers on using this method (Haiman & Loeb 1997, 
1998; Valageas & Silk 1999; Chiu & Ostriker 1999), but Madau et al. (1999) paper presents the 
clumping factor approach in its essence, in the simplest possible fiavor, and therefore is easy to 
reproduce. 



Figure |ll[ b now shows the most important quantity in the semi-analytical modeling: the 
porosity of the IGM Qhii as a function of time. The dashed line shows the approximate expression 
from equation (23) of Madau et al. (1999) paper, and the solid line refiects a more accurate 
calculation based on equation (21) of Madau et al. (1999) paper with the comoving emissivity 
extracted from the simulation. Again, the main conclusion of Madau et al. (1999) is confirmed, 
the overlap occurs when Qun ~ 1- 

I have just demonstrated, that my simulation agrees (at least on a semi-qualitative level) 
with semi-analytical models of Miralda-Escude et al. (1999) and Madau et al. (1999). Then how 
about the disagreement between the two? Clearly, if both agree with the same simulation, they 
also agree with each other. Let me know discuss the "differences" between the two semi-analytical 
approaches. 

First of all, if one looks carefully, one can notice that equation (23) of Madau et al. (1999) is 
indeed the ratio of e to i? using Miralda-Escude et al. (1999) notation, and the Miralda-Escude 
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et al. (1999) equation (3) e ^ R is indeed equation (24) of Madau et al. (1999). In other words, 
the dashed Une in Fig. [Tl] b is a ratio of the two curves in Fig. There is therefore no principal 
difference between the two approaches, but the advantage of Miralda-Escude et al. (1999) approach 
is that they propose a model for the clumping factor based on the realistic density distribution 
function, rather than make an ad hoc assumption about the time evolution of the clumping factor. 
The main disadvantage of their model is that it is not necessarily accurate quantitatively.^ 

There however still remains a question of the number of ionizing photons per baryon required 
to ionized the universe. It cannot be simultaneously one (Miralda-Escude et al. 1999) and many 
(Madau et al. 1999), can it? 

The irony of the situation is that it indeed can. The difference between these two approaches 
is mostly terminological, and it can be best understood if one considers the number of ionizing 
photons per baryon as a function of scale, N^/^iR) (this quantity cannot be rigorously defined 
mathematically, but has an intuitively clear physical sense). If Madau et al. (1999) count all 
photons emitted by all sources (in our case stars), i.e. they consider a quantity N^/ii{0) (here zero 
actually means the stellar surface), then Miralda-Escude et al. (1999) use N^/i){lGM) where the 
symbol IGM stands for the photons that escape local absorption, i.e. absorption in the gas bound 
to a given source, and corresponds to a range of scales of the order of 0.1 — 1 Mpc. The top panel 
of Fig. |l^ gives yet another quantity, N^/biRmin)^ where -Rj^i^ is the resolution of the simulation, 
and in my case it is 

^min = Y^^^^kpc 

in physical units. It is now easy to bridge a difference between Madau et al. (1999) and 
Miralda-Escude et al. (1999). The number A'^/j(0) is perhaps a factor of a few larger than 
^-y/biRmin) ~ 10 at 2 = 7, which gives a value for l/N^/ij{0) of the order of a few percent, in 
complete agreement with Madau et al. (1999). On the other hand, as has been discussed above, 
A^^/;,(IGM) should be considerably smaller than ^-y/biRmin)^ which leaves us with the conclusion 
of Miralda-Escude et al. (1999) that Ny/f){lGM) is of the order unity. 

There remains one final sticking point: a question of whether reionization is "fast" or "slow". 



But I have already mentioned in § p.l| that it is to a some degree a question of terminology as well. 
The whole process of reionization consists of three different stages: the "slow" (of the order of 
the Hubble time) pre-overlap, the "fast" (of the order of one tenth of the Hubble time) overlap, 
and the "slow" post-overlap. Miralda-Escude et al. (1999) focused on the post-overlap stage, but 
considered the whole process of reionization, and therefore concluded that reionization was slow. 
Madau et al. (1999) (and many other authors including Gnedin & Ostriker 1997) used the term 



^Miralda-Escude et al. (1999) model will work much better quantitatively for the case when the universe was 
reionized by a few bright quasars. In this case the amount of gas in the vicinity of sources is much smaller than in 
the case of stellar sources of reionization, when there are sources in a good fraction of all high density regions. 

^There may however still remain a difference, since Miralda-Escude et al. (1999) also claimed that the overlap is 
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"reionization" to label the overlap stage only, and thus concluded that reionization was fast. Thus, 
all the semi-analytical models give more-or-less similar predictions for the physics of reionization, 
and these predictions are reproduced by numerical simulations. 

3.5. Reionization and the Lyman-alpha forest 

While a question of reionization has its own astrophysical interest, reionization also leaves 
an imprint on the IGM at lower redshift, i.e. on the Lyman-alpha forest. One of the important 
properties of the forest is the so-called "effective equation of state" , i.e. a tight relationship between 
the gas density and temperature in the low density regime, which is usually can be approximated 
by a power-law over a limited range of densities around the cosmic mean, 

(Hui & Gnedin 1998). The time evolution of two parameters. To and 7, depends on the thermal 
history of the universe, and therefore is tightly coupled to the processes taking place during 
reionization. 

Usually, the evolution of the equation of state is calculated using the optically thin 
approximation and assuming that the radiation field is uniform (Hui & Gnedin 1998). This is 
clearly an approximation which needs to be verified. 

Figure [l^ shows the comparison of the "effective equation of state" from the simulation and 
the one computed in the optically thin approximation using the evolution of the spatially averaged 
radiation field extracted from the simulation (and thus both calculations have precisely the same 
mean J^, as a function of time). While the two calculations are quite different at high redshift, the 
optically thin approximation gives a reasonably accurate answer at z ^ 5. 

It is however important to notice, that the relationship between the gas density and 
temperature is more complicated at around the epoch of overlap and before. Figure |l^ shows 
the joint mass-weighted distributions of the gas density and temperature at four redshifts, 
corresponding to Fig. ^ and Fig. ^. One can see that the power-law fits to the "effective equation 
of state" give poor representation of the true temperature-density relation at 2; ^ 6. 

Another important consideration that severely affects our ability to model the Lyman-alpha 
forest is how homogeneous the ionizing background is. Figure ^ shows the mean and rms 
fluctuation in the local ionizing background as a function of gas density at three values of redshift. 
As can be seen, there is essentially no evolution in the spatial distribution of the ionizing radiation 



"gradual". It is not clear however without further elaboration of the Miralda-Escude et al. (1999) model whether 
this "gradual" overlap is actually "fast" enough to be compatible with other semi-analytical approaches and with 
this work (after all, the ionizing intensity has to increase by two to three orders of magnitude over a Hubble time or 
so), or whether there is a genuine disagreement between Miralda-Escude et al. (1999) and this work. 
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Fig. 12. — Evolution of the "effective equation of state" of the IGM, as represented by parameters To {bottom 
panel) and 7 {top panel), for the full numerical simulation {bold line) and the optically thin approximation 
{thin line). The shaded line marks the redshift above which the temperature-density correlation is not well 
established. 

in the low density regime, and gas with overdensities in excess of about 10 is subject to severe 
(more than 30%) inhomogeneities in the ionizing background. 

Since in the low density regime there exist a tight correlation between the column density of 
an absorption line and the gas density the absorption line originates in, it is possible to convert 
the gas density into the column density of the absorption line. Assuming J21 = 0.5 independently 
of redshift, and using equation (18) of Ricotti, Gnedin, Sz Shull (1999), I obtain: 

A^Hi ~ 1-13 X 10^^ cm-2(l + 5)^-^2 at z = 4, 

and 

Nm ~ 2.57 X 10^^ cm.~^{l + S)^-^^ at z = 3. 
Thus, if a 30% fluctuation in the ionizing background is still compatible with the background being 
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Fig. 13. — The joint mass-weighted distribution of the gas density and temperature in the production run 
N128X4_A at four different redshifts. Bold sohd lines show the power-law fits to the "effective equation of 
state" at the respective redshifts. 

uniform, only the Lyman-alpha forest with A?hi < 4 x lO-*^^ cm~^ at z = A and iYni < 1 x 10-^^ cm~^ 
at z = 3 can be modeled in a uniform ionizing background approximation. If an accuracy of 10% 
is required (which corresponds to 6 < 1), then these limits shrink further to 3 X lO^^cm"^ and 
8 X 10^^ cm~^ respectively. 

It is important to remember, that the production run simulation presented in this paper has 
still not converged numerically to below several tens of percent accuracy, and thus the numbers 
quoted above should be considered as only "within-a-factor-of-two" limits. Thus, it is likely to 
be safe to claim that theoretical models that assume a homogeneous ionizing background cannot 
reproduce the Lyman-alpha forest with the column density 

ArHi>10i5cm-2 
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Fig. 14. — Average {thin lines) and rms {bold lines) fluctuation in the ionizing background as a function of 
gas density at z = 5 {dotted lines), z — 4.5 {dashed lines), and z ~ A {solid lines). 

at a better than 30% level, and the forest with 

iVni > 10^^cm~2 

with a better than 10% accuracy for z between 3 and 4 on a line-by-line basis.0 



4. Conclusions 

Cosmological numerical simulations that incorporate the effects of radiative transfer show 
that the whole process of reionization can be separated into three main stages: on a "pre-overlap' 
stage, which occupies a Hubble time or so, the Hii regions expand into the neutral low density 



^Of course, these errors may average out for some of the global properties of the forest. 
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IGM, leaving high density protrusions behind. The high density regions in the vicinity of a source 
are being ionized as weU, albeit with a much slower rate. Pre-overlap culminates with the overlap, 
when the H li regions finally merge and the universe becomes transparent on a much larger scale 
than during the previous stage. The overlap epoch is characterized by a sharp rise (in about 10% 
of the Hubble time) in the level of the ionizing background and in the mean free path. At the 
"post-overlap" stage the universe consists of highly ionized low density gas, and still neutral high 
density gas, which is being gradually ionized as more and more ionizing photons become available. 
The boundary between the neutral and the ionized gas moves toward higher densities with time. 

In a currently fashionable CDM-type cosmological scenarios the whole process of reionization 
occupies a considerable part of the early evolution of the universe from z ~ 15 until z ~ 5. 

The aftermath of reionization has a profound effect on the IGM at lower redshift: the residual 
fluctuations in the ionizing background become quite significant at the column densities in the 
range from 10^^ cm~^ to 10^^ cm~^ (at z = 3 — 4) and may compromise theoretical models that do 
not take these fluctuations into account. 
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A. Local Optical Depth Approximation 

As was shown in Gnedin & Ostriker (1997, equation [B9]), the energy density of radiation per 
unit frequency Iu{x^,t) can be represented as 

Ux\t) = 7,(t)e^-'W---^(-^*) + / d3xi^^^lM^^l^^e--(-^-i'*), (Al) 

47rc J (x* — 

where I,y{t) is the volume average of I,^, Sy is the source function, Ti,{x'^,x\) is the optical depth 
between and x\, Ty{x'^,t) is the optical depth from a given point x* to an "average" point in 
the universe, obtained by solving the homogeneous radiative transfer equation with no source 
function, and fy{t) is the normalization constant obtained from the following condition: 

with brackets denoting the volume average. 

In order to implement a radiative transfer scheme into modern cosmological simulations, 
one needs to be able to solve equation ( [Al] ) in 0{Nb) operations (where Nb is the number of 
baryonic resolution elements, and I ignore a possible logarithmic multiplier, so that O^NBlogNs) 
still counts as 0{Nb))- However, equation (Al) in general includes two operations that are 
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more expensive than this count: (1) the one-point optical depth Tu{x^,t) requires 0{Nb x Nd) 
operations, where Nd is the number of directions, and in general is of order of n"^^ , and (2) 
direct evaluation of the integral over the source function requires 0{Nb x Ng) operations, where 
Ng is the number of sources. In the case under consideration, Ns is the number of stellar particles 
in the simulations, and is typically comparable to the number of baryonic resolution elements 
(perhaps, a factor of a few smaller, but it is the scaling that matters). Thus, direct (i.e. exact) 
implementation of equation (Al) is not feasible at the moment. ^ 



Thus, one has to use approximations to try to decrease the operation count in evaluating 



(Al). The "Local Optical Depth" approximation to evaluate the first term in (Al) was developed 
in Gnedin & Ostriker (1997) and is based on a simple observation that high precision in evaluating 
T^(a;*,t) is not required: if the optical depth is small, it does not matter much whether it is 10^^ 
or 10~^. If it is large, it again does not matter much whether it is 10^ or 10^. Only when r ~ 1 
it needs to be computed accurately, but the volume of space where this condition is achieved 
is very small, so at the end no large error is introduced even if r is computed with only an 
order-of-magnitude accuracy. Based on this consideration, the first ansatz is introduced in the 
solution of equation (|A1|): 

4"^n(")(x^)L(x*), (A2) 



a 



where index a runs over the list of species (in our case Hi, Hei, and Heii), and L is a characteristic 
length, which is taken to be the characteristic length of the density distribution, and I drop the 
time dependence (present in all terms) hereafter for the sake of simplified notation. Specifically, I 
adopt the following expression for L: 

L = I , ^ (A3) 

'a|Vlogp| + /3p |Alogp| + 7 |VlogXHi| 

where p is the baryon density, xhi is the neutral hydrogen fraction, a = 0.216 and (5 = 0.068 are 
constants |^ chosen to reproduce the correct result for the column density of a l/(r^ + r^) density 
distribution in two limits r ^ and r ^ oo and the last term is designed to work only for uniform 
density distribution, when L otherwise would be infinite. To this end I choose 

7 = 30 max(0, 1- p/p), (A4) 

where p is the average density of the universe and the factor 30 in front is chosen in order to satisfy 
the test D discussed above. Physically, this coefficient means that the approximation described in 
this paper spreads ionization fronts on average by a factor of \/30 ~ 5. 



*In the case of reionization by quasars, however, when the number of sources Ns is small enough, the integral in 



(Al) can indeed be computed directly in 0{Nb) operations (Abel et al 1999). 

^Note, that in Gnedin & Ostriker (1997) these constants are listed incorrectly because of a typo; they had proper 
values in the simulations described in that paper. 
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The Local Optical Depth approximation ensures that the first term in equation ( [All) is 
computed in 0{Nb) operations. 

Let me know turn to the second term. In general, integrals similar to the one in (Al) can only 
be computed in 0{Nb) operations with high quasi-Lagrangian resolution if they are convolutions. 
Thus, one of the ways to make this integral computable is to approximate it as a sum of several 
convolutions. For this purpose I divide all sources into two categories, "far" and "near" sources, 
and rewrite equation (|A1[) as 

Ux') = 7,e^''-^^("') + I^{x') + I^{x'). (A5) 

For a fluid element in a cosmological simulation, a source can be considered to be "far" if it sits 
in a separate clump. Then, the optical depth between the fluid element and the source can be 
approximated as a sum of one-point optical depths from two points: 

and, thus. 



This integral is a convolution, and can be computed in 0{Nb) operations by means of a standard 
P^M technique. More than that, since the SLH-P^M code already incorporates a P^M gravity 
solver, it has been straightforward to modify the existing solver to compute the (Al) term. 

The second, "near" term, //X (x*), can be considered to include everything that is left to 
include in equation ( [AS] ) after the approximation ( [A6D is incorporated. However, in order to 
overcome the 0{Nb x Ns) count, I must introduce another ansatz. For a "near" source, the 
two-point optical depth Tu{x'^-,x\) can be decomposed into Taylor series, 

Ty{x\x\) K, VjT„{x\x'){x^ - x{). 

However, even this is not enough since the factor ViT^{x''',x^) still prevents the integral in ( [Al[ ) 
from being a convolution. Thus, I introduce the second ansatz in the following form: 

jNr^i) = — f rf3^ /;(:"l)-5 e-^l"'-"il \l - e-^K)-^(-')l . (A7) 
47rc J (x* — x\)'^ L J 

Here 5 is a constant, independent of position (but it still may depend on time), and the integral 
now is a sum of two convolutions. The square bracket ensures that in the optically thin limit, 
r — > 0, the "near" term vanishes and the "far" term reduces to 



Attc J (x* — x\Y ' 

which is indeed the integral in equation (Al) in the limit r ^ 0. Thus, the proposed approximation 
becomes exact in the optically thin limit. 
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The fact that the factor 5 is the same for ah sources imphes that the ionization front speed 
is correct only "on average", and each particular ionization front may propagate with a speed 
different from the correct one, but they will merge at approximately the right moment. Of course, 
if there is only one ionization front in the simulation, its speed is computed correctly. 

From the argument presented above, it is clear that a good form for the factor S would be 
the following one: 

where thi{x^) is the one-point optical depth at the hydrogen ionization threshold, and the ratio 
thj I L represents the gradient of the optical depth with the desirable feature that it is independent 
of the poorly defined quantity L. 

The quantity C, which in part determines the speed of ionization fronts, needs to be fit to the 
tests. In other words, there is noting in the approximation ([A5- A8) which ensures photon number 



conservation, so the quantity C should be chosen so that this number is at least approximately 
conserved. 

Finally, a few words are in order about the frequency dependence. For stellar sources, 
considered in this paper, the shape of the source function is the same everywhere, so I can take it 
out of the integrals, 

Suix\t) = guPMsix\t), 

where is constant in space and time, and pMS is the mass density of massive stars. However, 
this still leaves the frequency dependence of the optical depth in and With the current 
computer capabilities it presents a considerable expense to perform the calculation of the integral 
in and at a sufficient number of frequency values. Therefore, I introduce the third ansatz 
in the described approximate scheme. I introduce three effective column densities N^f[\ where a 
again runs over Hi, Hei, and Heii, so that 

and analogously for l'^ . Then the integrals in and only need to be computed at zero 
frequency and at three threshold frequencies of Hi, Hei, and Heii. In the future, when computer 
power increases by another factor of 10 or so, it will be possible to get rid of the third ansatz and 
compute integrals in and over a large number of frequency bins. 



^°Tests show that uniform logarithmic sampUng of the frequency space should be at least 50 points per decade, 
which makes the total number of required frequency bins at least 200. 
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B. Testing the Approximation 

Obviously, the proposed approximation is so simple that it is highly unlikely that it will 
work in all circumstances. Thus, my goal is to make sure that it works in cosmologically relevant 
conditions. This can be achieved by appropriately "training" the scheme against spherically 



symmetric numerical solutions, which can be obtained from direct integration of equation (Al). 
For the test, I choose the following density distribution: 

p{r) = (Bl) 

This density distribution peaks at the center with p = pc, has a core radius of (taken to be 11 
comoving kiloparsecs in all cases), and falls off as r"" outside of the core radius. The specific 
form of the density law was chosen to make the coordinate transformation from real space to 
Lagrangian space g*, 



analytical. 



pd^x = pd^q, 



to \- 3 \ 3/{3-n) 



?,pcrl 



1/3 



This distribution is then embedded into the uniform mesh = at a point where the average 
density inside the sphere with radius is p, so that {p) = p over the computational box. Tests 
are then performed with the 64^ computational Lagrangian mesh and compared to the results of 
a spherically symmetric calculation. Spherically symmetric calculations are done with the same 
density profile and with 4000 shells between 0.1 kpc to 100 kpc. It has been verified that this 
number of shells is sufficient to achieve a 1% accuracy in the exact solution. 

As the first test, I adopt a n = 2 (p oc ar large radii) density distribution with pc/p = 10^ 
(maximum overdensity of 1000) embedded in the uniform radiation field with F = 10^^^ s^^ and 
F = lO-^'^s-i (which roughly correspond to the radiation field after and before the overlap) at 
z = 9. No point source is included in this test, so it designed to test only the expression for the 



optical depth (A2,A3). Figure 15 now shows the temperature and the ionization fraction for 
the exact and the approximate solutions taken at some arbitrary moments in time for the two 
values of the photoionization rate. As one can see, the neutral hydrogen profiles are reproduced 
reasonably well, except just after the ionization front, whereas temperature profiles are reproduced 
less accurate, with extra heating observed just after the ionization front. As I will show below, the 
same behavior is observed in the test cases with the source places at the center and the ionization 
front propagating outward. This feature has to be considered an error of the LOD approximation 
and cannot be easily repaired. 

Next, I perform a series of test with a source of ionizing radiation placed at the center of 
the spherically symmetric density distribution. The strength of the source is parametrized by the 
photoionization rate at the core radius Fc computed in the optically thin regime. 
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Fig. 15. — Distribution of the temperature (top panel) and neutral hydrogen fraction on linear (middle 
panel) and logarithmic (bottom panel) scale for a homogeneous radiation field test shining on the n = 2 
and pc/p = 10^ density distribution from outside at z = 9 with the photoionization rates of F = 10~^^s~^ 
and F = 10^^'* s^^ taken at an arbitrary moments in time. Solid lines show the exact spherically symmetric 
solutions, and symbols mark the approximate 3D solutions, with solid squares corresponding to F = 10^^^ s^^ 
and open squares corresponding to F = 10^^"* s~^. 

The seven main tests performed are listed in Table The first five tests include the density 
distribution that falls off as at large radii, the test F verifies the propagation of the ionization 
front over the uniform density field {pc/p = 1); and the last test G includes the density law as 
the intermediate case between the tests A-E and F. That test is not used to train the approximate 
scheme, but rather to verify the performance of the approximation in a substantially different 
arrangement. 

Additional tests have been performed to verify numerical convergence and other technical 
issues such as dependence on the cell size. 

As the results of fitting the approximate scheme to the exact solutions of the first six tests. 
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Table 2. Parameters of Test Cases 



Case 


z 


Tea' (s-i) 


Pc/p 


n 


*o (yr) 


A 


4 


10-14 


10^ 


2 


1 X 10^ 


B 


4 


10-12 


10^ 


2 


2 X 10"^ 


C 


4 


10-10 


10^ 


2 


4 X 102 


D 


9 


10-14 


10^ 


2 


1 X 10^ 


E 


9 


10-12 


10^ 


2 


2 X 10^ 


F 


9 


10-14 


1 




1 X 10^ 


G 


9 


10-12 


40 


1 


1.5 X lO'' 



the following expression for the coefficient C is derived: 



where 



and 



r = max [O, (r - ^R(T)ne)j 



where T is the photoionization rate, R is the recombination coefficient, and is the electron 
number density. The quantity C thus depends on the time integral of the volume averaged 
photoionization rate (with a correction for recombinations), and the factor a^^^ in the denominator 
accounts for the redshift of ionizing photons. 



Equation (|B2|) thus completes the approximation. I would like to stress here that equation 
(B2) is a purely empirical fit and may have no physical basis behind it. It appears physically 
plausible since it is based on the number of ionizing photons emitted over the course of action, but 
its specific form cannot be "deduced" from the first principles (in large part due to the fact that 
the Local Optical Depth approximation is based on two ansatzes). 

We can now investigate how the developed approximation works in test cases. Figure ^ 
shows the comparison between the exact and approximate calculations of the evolution of the 
ionization front in all seven test cases. The agreement is good everywhere except cases D and E at 
early times, but this is due to a finite resolution of a numerical simulation. Deviations of the order 
of 20% are also observed in case A at late times. The case G shows in general a worse agreement, 
because this test was not used to fit the approximate scheme to the exact solutions. It can thus be 
used to demonstrate the level of accuracy of the LOD approximation for arrangements different 
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Fig. 16. — Evolution of the ionization front in seven test cases. The upper panel shows cases A (dotted 
line), B (short-dashed line), and C [long-dashed line). The lower panel shows cases cases D (dotted line), E 
(short-dashed line), F (long-dashed line), and G (dot-dashed line). Thin lines mark the exact solution, and 
bold lines show the Local Optical Depth approximation. Time axis is scaled arbitrarily, with the scaling 
factor listed in Table ^j. The thin solid line shows the Shapiro analytical solution for the test F (homogeneous 
density field). 



from those used to "train" the approximation. The solid thin line in Fig. 16 shows the analytical 
approximation of Shapiro (1986) and Shapiro & Giroux (1987) for the case F (the homogeneous 
density field). The small deviation observed at early times is due to the fact that the thickness of 
the ionization front in this regime is not negligible compared to the size of the Hii region, and the 
Shapiro solution becomes invalid. 

The question of numerical resolution in any 3D treatment of the radiative transfer is of utmost 
importance. Gravitational collapse proceeds from low density to high density, and therefore even 
with poor resolution low density structures can be reproduced reliably. The situation is just the 
opposite with the radiative transfer: the ionization front moves from high density regions into 
low density regions. Thus, during the first time-step the ionization front has to be fully resolved, 
or it will never leave the resolution element it originated in (unless it is specifically followed on 
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Fig. 17a 



Fig. 17b 



Fig. 17. — Distribution of the temperature {top panel) and neutral hydrogen fraction on hnear {middle 
panel) and logaritlimic {bottom panel) scale at five different times (as labeled by different lines) in the exact 
solution {thin lines) and the approximate solution {thick lines) for cases A (panel a), D (panel b), and B 
(panel c). 

inside that resolution element). Thus, to avoid the problem of ionization front getting "stuck" due 
to lack of resolution, I impose the condition that Se > 2, where e is the gravitational softening 
length. In test cases E and F this leads to the ionization front deviating from the exact solution 
at earlier times, but catching up later, when it becomes fully resolved in the approximate scheme. 



Figure 17 shows in three panels time-evolution of the temperature and the neutral hydrogen 
fraction for cases A, B, and D. Since the ionization front in the approximate solution is not sharp, 
I choose the location where hydrogen is 50% ionized as a location of the ionization front in the 



approximate solution for the purpose of producing Fig. 16 



Several defects of the approximate scheme are apparent in Fig. First, the ionization front 
is spread over a considerable distance, but for cosmological purposes this distance is still much 
smaller than the correlation length, and thus this does not make the big impact on the dynamics of 



the gas. Second, at late times in Fig. 17a the approximate scheme significantly overestimates the 
degree of ionization, due to the fact that the optical depth between the position of the ionization 
front and the source is of the order of unity, and here the Local Optical Depth approximation 



is bound to fail. Third, in Fig. 17a the temperature ahead of the front is higher than it should 



be, whereas in Fig. 17b it is lower. This is due to the fact that leaking of high-energy photons 
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Fig. 18 

Fig. 18. — Distribution of the temperature {top panel) and neutral hydrogen fraction on hnear {middle 
panel) and logarithmic {bottom panel) scale in test B in the exact solution {thin lines) and the approximate 
solution {thick lines). Dotted lines show the full solution, and dashed lines show the solution with cooling 
disabled. Cooling along the ionization front leads to an incorrect prediction for the post-front temperature. 

across the front is not captured correctly by the approximation. Finally, the most severe defect is 
apparent in Fig. |l7| g: the temperature inside the front in the approximate solution is some 50% 
lower than it should be. This is due to the fact that the ionization front is not sharp, and the gas 
is able to cool while being ionized, whereas in the exact solution the neutral fraction drops from 1 
to 10~^ so quickly, that no appreciable cooling occurs. To illustrate this further I show in Figure 
18 the first line from Fig. ^b together with the exact and approximate solution for the same test 
case but with disabled cooling. Without cooling, the approximate solution predicts the post-front 
temperature right, but with cooling it allows a considerable loss of energy. 

This defect however is not specific to the Local Optical Depth approximation, but obviously 
will be present in any scheme that does not resolve the ionization front. It is possible to modify 
the SLH technique in such a way as to allow the mesh to deform appropriately and to achieve 
higher resolution at the ionization front, but this work is well beyond the scope of this paper. 

Another significant defect of the Local Optical Depth approximation is the absence of 
shadowing. Since LOD approximation does not include proper ray tracing, the ionizing flux from a 
given source in a given fluid element is only influenced by the opacity in the vicinity of the source 
and the fluid element, and therefore a shadowing effect of a self-shielded neutral cloud on the line 
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of sight from the fluid element to the somxe would not be reproduced. Fortunately, shadowing 
is not important during cosmological reionization: a self-shielded cloud with the size of lOkpc 
moving with the velocity of 100 km/ s produces a shadowing effect only during a time interval of 
100 Myr, which is more than two orders of magnitude shorter than the recombination time in the 
low density IGM. Thus, the ionization and the thermal state of the IGM would not be affected by 
the passage of a cloud shadow. 
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